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This review concentrates on the two principle methods used to evolve nuclear abun- 
qv | dances within astrophysical simulations, evolution via rate equations and via equi- 

~~^ libria. Because in general the rate equations in nucleosynthetic applications form 

an extraordinarily stiff system, implicit methods have proven mandatory, leading 
1. ', to the need to solve moderately sized matrix equations. Efforts to improve the per- 

formance of such rate equation methods are focused on efficient solution of these 
matrix equations, by making best use of the sparseness of these matrices. Recent 
work to produce hybrid schemes which use local equilibria to reduce the compu- 
tational cost of the rate equations is also discussed. Such schemes offer significant 
^ improvements in the speed of reaction networks and are accurate under circum- 

3 . stances where calculations with complete equilibrium fail. 



1 Introduction 



By the second half of the last century, work by Helmholtz, Kelvin and others 
made it clear that neither gravity nor any other then known energy source 
could account for the age of the sun and solar system, as determined by ge- 
ological measurements. Following quickly on the heels of Rutherford's 1911 
discovery of the atomic nucleus, Eddington and others suggested that nuclear 
transmutations might be the remedy to this quandary. With the burgeon- 
ing knowledge of the properties of nuclei and nuclear reactions in the 1930s, 
1940s and 1950s, came a growing understanding of the role that individual 
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Fig. 1. The abundances of isotopes in the solar system as a function of atomic mass 
[3,4]. The abundances are normalized so that the total abundance of silicon is 10 6 . 

nuclear reaction played in the synthesis of the elements. In 1957, Burbidge, 
Burbidge, Fowler & Hoyle [1] and Cameron [2] wove these threads into a co- 
hesive theory of nucleosynthesis, demonstrating how the solar isotopic abun- 
dances (displayed in Fig. 1) bore the fingerprints of their astrophysical origins. 
Today, investigations refine our answers to these same two questions, how are 
the elements that make up our universe formed, and how do these nuclear 
transformations, and the energy they release, affect their astrophysical hosts. 

In this article, we will concentrate on summarizing the two basic numerical 
methods used in nucleosynthesis studies, the tracking of nuclear transmuta- 
tions via rate equations and via equilibria. We will also briefly discuss work 
which seeks to meld these methods together in order to overcome the limita- 
tions of each. To properly orient readers unfamiliar with nuclear astrophysics 
and to briefly describe the differing physical conditions which influence the 
optimal choice of abundance evolution method, we begin with a brief intro- 
duction to the background astrophysics (§2), before discussing the form that 
the rate equations take (§3). In §4 we will discuss the difficulties inherent in 
solving these rate equations. §5 describes the equations of nuclear equilibria 
as well as the limitations of their use. Finally in §6 we will discuss hybrid 
schemes which seek to use local equilibria to simplify the rate equations. 



2 Nuclear Processes in Astrophysics 



In general, nucleosynthesis calculations divide into two categories, (1) nu- 
cleosynthesis during the hydrostatic burning stages of stellar evolution and 
(2) nucleosynthesis in explosive events. The critical distinction between the 
categories is the timescale over which transmutations occur. In hydrostatic 
burning, the release of nuclear energy occurs at a rate balancing the loss of 
energy via radiation and neutrinos (the photon and neutrino luminosities), 
providing continuous pressure support against the star's self-gravity This re- 
sults in slow burning at relatively low temperatures and densities. In explosive 
nucleosynthesis the timescale and thermodynamic conditions are determined 
by the hydrodynamics. For example, in many cases of interest, a detonation 
shock heats the material and expels it outward, to cool adiabatically as it 
expands. In such a case, the limiting timescale is that of the hydrodynamic 
expansion, not those intrinsic to the nuclear processes. Though constrained 
from an in depth discussion by our concentration on the computational meth- 
ods, in this section we will outline the physics of the many burning stages 
and in particular the way that the physics influences the choice of method for 
computing the abundance changes. 



2. 1 Hydrostatic Burning Stages in Stellar Evolution 



A star shines bright because, deep in its core, energy released by thermonu- 
clear reactions balances the star's self-gravity Table 1 lists the series of burning 
stages which make up a star's life, with their representative conditions for a 
star like the sun and a star 20 times more massive (M=20 M & ). With the con- 
sumption in turn of each nuclear fuel, the inexorable squeeze of gravity creates 
higher temperatures and densities. The initiation of each subsequent stage re- 
quires these progressively higher temperatures and densities to overcome the 
increasing Coulomb repulsion of the reactants. Around the core, layers of pro- 
gressively lighter matter echo prior central burning stages. A star's mass is 
the most important determinant of its final destiny, with less massive stars 
not progressing through the full sequence of burning stages. For stars with 
masses less than 8 M , helium burning is the final burning stage, because 
the packing of electrons into a degenerate Fermi-Dirac configuration provides 
sufficient pressure support to prevent further contraction. Instead, the remain- 
ing envelope is driven off, forming a planetary nebula, and leaving the bare 
core. Such cooling bare cores, composed of C & O in this case, or He or O, 
Ne & Mg in others, are termed white dwarf stars. As examination of Table 1 
reveals, in less massive stars the individual burning stages which occur do so 
at higher density, lower temperature and over much longer timescales. For a 
more exhaustive description of stellar evolution consult, e.g., [5-7] 



Table 1 

Stellar Burning Stages 
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Hydrogen 
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PP chain 
Triple a 


For a 20 M star [9] 



Hydrogen 5.6 

Helium 9.4 xlO 2 

Carbon 2.7xl0 5 

Neon 4.0 xlO 6 

Oxygen 6.0xl0 6 

Silicon 4.9 xlO 7 



0.040 

0.19 

0.81 

1.7 

2.1 

3.7 



l.OxlO 7 2.7xl0 38 - 

9.5xl0 5 5.3xl0 38 < xlO 36 

3.0xl0 2 4.3xl0 38 7.4xl0 39 

0.4 4.4xl0 38 1.2xl0 43 

0.5 4.4xl0 38 7.4xl0 43 

0.01 4.4xl0 38 3.1 xlO 45 



CNO Cycle 
Triple a 

12 C + 12 C _ 20 Ne + a 

20 Ne + 7 -> 16 + a 
i6 + i6 ^28 Si + a 

28 Si + 7a -» 56 Ni 



Because the hydrostatic burning timescales are long compared to beta-decay 
half-lives (with a few exceptions of long-lived unstable nuclei), nuclei can de- 
cay back to stability before undergoing the next reaction. As a result, most 
reactions in hydrostatic burning stages proceed through stable nuclei. As in- 
dicated in Table 1, stars transmute hydrogen into helium via two alternative 
reaction sequences; the PP-chains, initiated by the conversion of two protons 
to a positron, a neutrino and a deuteron, (commonly notated 1 H(p, e + z/) 2 H), 
and the CNO cycle, which converts 1 H into 4 He by a sequence of catalytic 
(p, 7) and (p,a) reactions on pre-existing C, N, and O nuclei and subsequent 
beta-decays. He-burning follows H-burning, converting 4 He into 12 C via the 
triple-alpha process, which we will discuss further in §6. A portion (larger in 
more massive stars) of these 12 C nuclei capture an additional a-particle to 
form 16 0. With the exhaustion of 4 He, contraction continues until conditions 
are sufficient to fuse pairs of 12 C nuclei, producing 20 Ne and a small fraction of 
23 Na. Depletion of 12 C is followed by further contraction. Before temperatures 
become sufficient for 16 nuclei to fuse, the thermal photon bath becomes en- 
ergetic enough to photodisintegrate 20 Ne (see §3.1 and Eq. 8 for more details), 
freeing an a-particle which can be captured by other 20 Ne nuclei forming 24 Mg. 
Once 20 Ne is exhausted, continued contraction raises the temperature until it 
is sufficient for 16 to fuse, producing 28 Si and 31 P. 

Si-burning, like Ne-burning, is initiated by photodisintegration reactions which 
then provide the light particles needed for capture reactions. Temperatures 
sufficient to photodisintegrate Si are also sufficient to photodisintegrate all 
other nuclei and permit charged particle captures (and, in any case, neutron 
captures, where there is no Coulomb barrier to overcome). This leaves many 



individual reactions in a chemical equilibrium where reactions are balanced by 
their inverses. If each of the important reactions connecting two species is in 
chemical equilibrium, then the relative abundances of these species will obey 
an equilibrium distribution. As the most bound nuclei, and therefore most 
resistant to photodisintegration, isotopes of Fe and Ni dominate the products 
of silicon burning, producing the Fe-peak seen near A = 60 in Fig 1. We 
will discuss this equilibrium further in §5. For a more complete description 
of silicon burning, see [10,11]. It is important to note that because of the 
large photodisintegration fluxes which nearly balance their reverse capture 
reactions, the effective rate at which silicon is burned is as much as 10 5 times 
slower than the individual reaction rates would indicate. This near balance of 
reaction rates during silicon burning will also prove important in our discussion 
of numerical methods in §6. More extensive overviews of the major and minor 
reaction sequences in all burning stages from helium to silicon burning in 
massive stars is given in [7,12-15]. 

An additional hydrostatic nucleosynthetic process is the s(low neutron capture)- 
process. While energetically unimportant, this slow neutron capture process 
leads to the build-up of the (small) abundances of roughly half of the heavy 
elements. During core and shell He-burning, (a,n) reactions on neutron rich 
nuclei provide a source of free neutrons, initiating a series of neutron captures 
and /9-decays, starting on pre-existing medium and heavy nuclei (up to Fe), 
synthesizing nuclei up to Pb and Bi. Such reactions, occurring under tem- 
perature conditions where photodisintegration are unimportant, approach a 
steady flow, wherein the flux of A_1 Z(n, 7) A Z equals the flux of A Z(n, 7) A+1 Z, 
as long as the relevant reaction timescale is short in comparison to the process 
timescale. For general overviews of the s-process see [16-19]. From a numer- 
ical perspective, in this and several other cases, it is often advantageous to 
couple only the dominant energy producing reactions to the hydrodynamics 
and perform detailed nucleosynthesis, which can require a very large nuclear 
network, in a post-processing fashion. 



2.2 Explosive Burning 



While hydrostatic sources are capable of producing many of the isotopes shown 
in Fig. 1, most of these products are trapped deep in the potential well of their 
parent star or white dwarf. Explosions are necessary to liberate the transmuted 
nuclei from this gravitational embrace. In massive stars, the formation of the 
iron core during silicon burning marks the end of nuclear energy generation 
in the core, as nuclei more massive than the iron peak nuclei are less bound. 
Ultimately, gravitational contraction turns into collapse, which is halted sud- 
denly by degenerate nucleon pressure when the matter density approaches or 
exceeds that of the atomic nucleus (~ 10 14 gcm~ 3 ). Infalling material from 



overlying layers bounce off of this newly formed proto-neutron star, sending 
a shockwave outward. Though this shock soon stalls, it is re-invigorated by 
neutrinos carrying off the gravitational energy released in the formation of 
the proto-neutron star (see, e.g., [20, this volume]). This re-energized shock 
propagates out of the star, unbinding much of the overlying layers and driv- 
ing them into space with velocities of thousands of km/s, and producing a 
(core-collapse) supernova [9,21,22]. This passing shock also causes further nu- 
cleosynthesis to occur as it passes through the layered composition of the star, 
raising temperatures to several GK. Many of the hydrostatic burning processes 
discussed in §2.1 occur under these explosive conditions but at much higher 
temperatures and over much shorter timescales. With little regard to the initial 
composition, the burning stage which determines the nucleosynthetic outcome 
is that whose timescale at the peak temperature is comparable to the hydro- 
dynamic timescale (see, e.g., [23-25]). Much more detail on the products of 
explosive burning in core collapse supernova is available in the literature (see, 
e.g., [7,14,15,25-27]). Because these explosive burning stages are responsible 
for producing many nuclei with masses between 16-70, detailed modeling is 
important. However such modeling is complicated by the strongly varying hy- 
drodynamic conditions. The computational feasibility of performing detailed 
nucleosynthesis calculations (particularly explosive silicon burning) within the 
latest generation of multi-dimensional hydrodynamic models is greatly aided 
by techniques like those we will discuss in §6.1 

The explosive analogue of the s-process is the r(apid neutron capture) -process 
which requires larger neutron concentrations. Conditions suitable for r-process 
nucleosynthesis can occur in the decompression of neutron star matter [28-32] 
or in the innermost regions of the core-collapse supernova ejecta [33,34]. As 
a result of the electron and neutrino captures this material has experienced, 
there are more neutrons than protons. With equal numbers of protons and 
neutrons tied up in 4 He, the remainder provide the required neutron to heavy 
seed nucleus ratio. Calculations independent of the specific astrophysical site 
[35,36], performed with the goal of reproducing the solar abundance pattern 
of heavy elements, show that extremely unstable nuclei close to the neutron 
drip line are produced and beta-decay timescales can be short in comparison 
to the process timescales. Because of the large number of nuclei involved, r- 
process calculations are among the most numerically expensive, however as 
we will discuss in §4 & §6, partial equilibria and/or the slow variation of the 
light particle abundances can be employed to simplify the calculation. 

A second type of supernova is the thermonuclear supernova, which occurs 
when explosive carbon burning is ignited in the center of an accreting white 
dwarf, either as a result of the white dwarf exceeding the maximum mass 
(~ IAMq) which can be supported by electron degeneracy pressure (see, e.g., 
[37-41]) or due to compression resulting from the detonation of an accreted 
He layer (see, e.g., [42,43]). In either case, a flame front propagates outward 



disrupting the white dwarf and leaving a composition dominated by iron peak 
and intermediate mass nuclei. Computationally, in addition to sharing the 
complications discussed for explosive burning in core collapse supernovae, ex- 
plosive nucleosynthesis in these thermonuclear supernovae is also the source 
of the energy which powers the explosion. Thus accurate hydrodynamic sim- 
ulation requires at least the inclusion of an accurate means to calculate the 
rate of thermonuclear energy release. 

If a white dwarf accretes hydrogen (typically via mass transfer from a binary 
companion) slowly enough, a layer of unburned hydrogen will build up on 
the surface of the white dwarf. Once the density of this layer is sufficient 
to ignite the hydrogen, a nova results (see, e.g., [44-46]). The degeneracy of 
the material and the steep gravitational gradients at the surface of the white 
dwarf, result in explosive hydrogen burning via the "hot" or [3 -limited CNO 
cycle [47,48], releasing 10 46 — 10 47 ergs over timescales of 100-1000 s, with 
peak temperatures reaching 0.2 — 0.3 GK. A neutron star can also accrete 
mass from a companion in a similar fashion, building up a layer of hydrogen 
which explodes to produce an X-ray burst (see, e.g., [49,50]). Because the 
neutron stars surface gravitational gradients are even stronger, the timescales 
are shorter (1-10 s) and even higher peak temperatures (1 — 2 GK) are reached, 
enabling proton capture up to and beyond the Fe-peak, the r(apid) p(roton 
capture) -process [47,51]. However the size of the hydrogen layer is much smaller 
(as small as 1O~ 12 M ) and, as a result, so is the energy release, which is 
typically 10 39 — 10 40 ergs. The thermonuclear source of the explosive energy 
in novae and X-ray burst, as well as the similarity in the convective and 
nuclear timescales, require that hydrodynamic simulations of these objects 
include large networks. Fortunately, approximations (which we will discuss in 
§6) exist which can greatly reduce the computational cost. 

An additional form of explosive burning occurred as the universe expanded 
from its primordial Big Bang. Since its origin, the universe has been expand- 
ing, cooling from the initially extreme temperatures. At the earliest times, 
the populations of all kinds of sub-atomic particles were equilibrated. Even- 
tually, continued cooling and the freezeout of this equilibrium left only a few 
neutrons and protons (approximately one baryon per billion photons) with 
approximately 7 protons per neutron. These nucleons largely remained free 
because the small Q-value of the reaction n + p v^ 2 H + 7 permits chemical 
equilibrium to persist to temperatures of 1 GK, keeping the abundance of deu- 
terium very small. By the time the universe has cooled to 1 GK, the expansion 
has reduced the density to ~ 10~ 5 gcm -3 , a density small enough that, even 
with an expansion timescale of days, the resulting tiny flux through the 3a 
reaction did not produce significant amounts of heavy elements. Instead Big 
Bang nucleosynthesis is responsible for the production of the lightest elements; 
1 H, 4 He, 2 H, 3 He and 7 Li [52,53]. In principle, the small number of affected 
species and the simple hydrodynamic evolution make Big Bang nucleosyn- 



thesis the least computationally challenging. However the greater accuracy of 
the relevant nuclear reactions, as well as the important limits that Big Bang 
nucleosynthesis calculations place on cosmologically important factors, have 
resulted in a much stronger drive to high precision in Big Bang nucleosynthesis 
calculations [54,55]. 



3 Thermonuclear Reaction Networks 



From the discussion in the preceding section, it is clear that nuclear abun- 
dances in many cases obey equilibrium distributions. But the general case 
requires the evolution of nuclear abundances via a nuclear reaction network. 
Composed of a system of first order differential equations, the nuclear reac- 
tion network has sink and source terms representing each of the many nuclear 
reactions involved. Prior to discussing the numerical difficulties posed by the 
nuclear network, it is necessary to understand the sets of equations we are at- 
tempting to solve. To this end, we present a brief overview of the thermonuclear 
reaction rates of interest and how these rates are assembled into the differ- 
ential equations we must ultimately solve. For more detailed information, we 
refer the reader to a number of more complete discussions [5,7,56-58] which 
can be found in the literature. We will end this section by briefly discussing 
the coupling of nucleosynthesis with hydrodynamics. 



3. 1 Thermonuclear Reaction Rates 



There are a large number of types of nuclear reactions which are of astrophys- 
ical interest. In addition to the emission or absorption of nuclei and nucleons, 
nuclear reactions can involve the emission or absorption of photons (7 -rays) 
and leptons (electrons, neutrinos, and their anti-particles). As a result, nu- 
clear reactions involve three of the four fundamental forces, the nuclear strong, 
electromagnetic and nuclear weak forces. Reactions involving leptons (termed 
weak interactions) proceed much more slowly than those involving only nu- 
cleons and photons. However these reactions are still important, as only weak 
interactions can change the global ratio of protons to neutrons. 

The most basic piece of information about any nuclear reaction is the nuclear 
cross section. The cross section for a reaction between target j and projectile 
k is defined by 

number of reactions target _1 sec _1 r/rij . . 

a = '" ' ' ' 



flux of incoming projectiles n^v 



The second equality holds when the relative velocity between targets of num- 
ber density rij and projectiles of number density n k is constant and has the 
value v. Then r, the number of reactions per cm 3 and sec, can be expressed 
as r = avrijUk- More generally, the targets and projectiles have distributions 
of velocities, in which case r is given by 

r j,k = / cr(\vj - v k \)\vj - v k \d 3 njd 3 n k . (2) 



The evaluation of this integral depends on the types of particles and distri- 
butions which are involved. For nuclei j and k in an astrophysical plasma, 
Maxwell-Boltzmann statistics generally apply, thus 

rf 3 n = n( ^_)3/2 ( _^l )d 3 (3) 

K 27ik B T J PV 2k B T J ' v ; 



allowing rij and n k to be moved outside of the integral. Eq. 2 can then be writ- 
ten as r^fc = (cv)j k rijnk, where (av) is the velocity integrated cross section. 
Equivalently, one can express the reaction rate in terms of a mean lifetime of 
particle j against destruction by particle k, 

r k {j) = -j-^ (4) 



For thermonuclear reactions, these integrated cross sections have the form 

[5,59] 

oo 

(j,k) = (av) jk = (_)V2 (fcBT) -s/2 [ Ea(E)eM-E/k B T)dE, (5) 

■" fin J 



where \x denotes the reduced mass of the target-projectile system, E the center 
of mass energy, T the temperature and k B is Boltzmann's constant. 

Experimental measurements and theoretical predictions for these reaction 
rates provide the data input necessary for nuclear networks. While detailed 
discussion of individual rates is beyond the scope of this article, the interested 
reader is directed to the following reviews. Experimental nuclear rates have 
been reviewed in detail by [56,60,61]. The most recent experimental charged 
particle rate compilations are the ones by [62,63]. Experimental neutron cap- 
ture cross sections are summarized by [64,65,18]. Rates for unstable (light) nu- 
clei are given, for example, by [51,53,66-69]. For the vast number of medium 
and heavy nuclei which exhibit a high density of excited states at capture 
energies, Hauser-Feshbach (statistical model) calculations are applicable. The 



most recent compilations were provided by [70-73]. Improvements in level den- 
sities [74] , alpha potentials, and the consistent treatment of isospin mixing will 
lead to the next generation of theoretical rate predictions [75,76]. 

In practice, these experimental and theoretical reaction rates are determined 
for bare nuclei, while in astrophysical plasmas, these reactions occur among a 
background of other nuclei and electrons. As a result of this background, the 
reacting nuclei experience a Coulomb repulsion modified from that of bare nu- 
clei. For high densities and/or low temperatures, the effects of this screening of 
reactions becomes very important. Under most conditions (with no n- vanishing 
temperatures) the generalized reaction rate integral can be separated into the 
traditional expression without screening [Eq. 5] and a screening factor, 

(J, k)* = fscr(Z j: Z k , p, T, 7ii)(j, k). (6) 



This screening factor is dependent on the charge of the involved particles, the 
density, temperature, and the composition of the plasma. For more details 
on the form of f scr , see, e.g., [77-81]. At high densities and low temperatures 
screening factors can enhance reactions by many orders of magnitude and 
lead to pycnonuclear ignition. In the extreme case of very low temperatures, 
where reactions are only possible via ground state oscillations of the nuclei 
in a Coulomb lattice, Eq. 6 breaks down, because it was derived under the 
assumption of a Boltzmann distribution (for recent references, see [81,82]). 

When particle k in Eq. 2 is a photon, the distribution d 3 n k is given by the 
Plank distribution, 

^ = ^e W (E/k B T)-l dE - ■ (7) 

Furthermore, the relative velocity is always c and thus the integral is separable, 
simplifying to 

fi — c ttz / ,,-, r, — ^ dE~ = A,--,(T)n,-. (8) 

3 n 2 (cTif J exp(£ 7 /A; B T) - 1 ' jM ! 3 K ' 



In practice there is, however, no need to directly evaluate the photodisinte- 
gration cross sections, because they can be expressed by detailed balance in 
terms of the capture cross sections for the inverse reaction, l + m — > j + 7 



Ai, 7 (T) = (^)(^) 3/2 ( ! ^) 3/2 (/,m)exp(-g, m /A; B T). (9) 
This expression depends on the partition functions, Gk = S?(2 J%+1) exp(— Ei/ksT) 
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(which account for the populations of the excited states of the nucleus), the 
mass numbers, A, the temperature T, the inverse reaction rate (l,m), and the 
reaction Q- value (the energy released by the reaction), Q im = (mi + m m — 
mj)c 2 . Since photodisintegrations are endoergic, their rates are vanishingly 
small until sufficient photons exist in the high energy tail of the Planck distri- 
bution with energies > Qi m . As a rule of thumb this requires TmQi m /30kB- 

A procedure similar to that for Eq. 8 applies to captures of electrons by nuclei. 
Because the electron is 1836 times less massive than a nucleon, the velocity of 
the nucleus j in the center of mass system is negligible in comparison to the 
electron velocity (\vj — v e \ ~ \v e \). In the neutral, completely ionized plasmas 
typical of the astrophysical sites of nucleosynthesis, the electron number den- 
sity, n e , is equal to the total density of protons in nuclei, J2i ZiUi. However in 
many of these astrophysical settings the electrons are at least partially degen- 
erate, therefore the electron distribution cannot be assumed to be Maxwellian. 
Instead the capture cross section has to be integrated over a Boltzmann, par- 
tially degenerate, or degenerate Fermi distribution of electrons, depending on 
the astrophysical conditions. The resulting electron capture rates are functions 
of T and n e , 

Vj = \ jt e(T,n e )nj. (10) 

Similar equations apply for the capture of positrons which are in thermal 
equilibrium with photons, electrons, and nuclei. Electron and positron capture 
calculations have been performed by [83-85] for a large variety of nuclei with 
mass numbers between A=20 and A=60. For improvements and application 
to heavier nuclei see also [86-89]. 

For normal decays, like beta or alpha decays, with a characteristic half-life 7*1/2, 
Eq. 8 & 10 also applies, with the decay constant Xj = ln2/ri/2- In addition 
to innumerable experimental half-life determinations, beta-decay half-lives for 
unstable nuclei have been predicted by [90,91] & [92, including temperature 
effects] . More recently, estimates have been made with improved quasi particle 
RPA calculations [93-100]. 

At high densities (p ~ 10 13 gcm~ 3 ), even though the size of the neutrino 
scattering cross section on nuclei and electrons is very small, enough scattering 
events occur to thermalize the neutrino distribution. Under such conditions the 
inverse process to electron capture (neutrino capture) can occur in significant 
numbers and the neutrino capture rate can be expressed in a form similar to 
Eqs. 8 & 10 by integrating over the thermal neutrino distribution (e.g. [101]). 
Inelastic neutrino scattering on nuclei can also be expressed in this form. 
The latter can cause particle emission, similar to photodisintegration (e.g. 
[102-106]). The calculation of these rates can be further complicated by the 
neutrinos not being in thermal equilibrium with the local environment. When 
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thermal equilibrium among neutrinos was established at a different location, 
then the neutrino distribution might be characterized by a chemical potential 
and a temperature different from the local values. Otherwise, the neutrino 
distribution must be evolved in detail [20, this volume]. 



3.2 Thermonuclear Rate Equations 



The large number of reaction types discussed in §3.1 can be divided into 3 func- 
tional categories based on the number of reactants which are nuclei. The re- 
actions involving a single nucleus, which include decays, electron and positron 
captures, photodisintegrations, and neutrino induced reactions, depend on the 
number density of only the target species. For reaction involving two nuclei, 
the reaction rate depends on the number densities of both target and pro- 
jectile nuclei. There are also a few important three-particle process, like the 
triple-a process discussed in §6, which are commonly successive captures with 
an intermediate unstable target (see, e.g., [107,108]). Using an equilibrium 
abundance for the unstable intermediate, the contributions of these reactions 
are commonly written in the form of a three-particle processes, depending 
on a trio of number densities. Grouping reactions by these 3 functional cate- 
gories, the time derivatives of the number densities of each nuclear species in 
an astrophysical plasma can be written in terms of the reaction rates, r, as 



dt 



= E^Vh + EK^, k + £A/j M r jVM , (ii) 

p=const j j,k j,k,l 



where the three sums are over reactions which produce or destroy a nucleus 
of species % with 1, 2 & 3 reactant nuclei, respectively. The M s provide for 
proper accounting of numbers of nuclei and are given by: M] = iVj, N]j k = 
Ni/UT=i \N jm \l, and H}^ = N i /U^ =1 \N jm \\. The N' t s can be positive or 
negative numbers that specify how many particles of species % are created or 
destroyed in a reaction, while the denominators, including factorials, run over 
the n m different species destroyed in the reaction and avoid double counting 
of the number of reactions when identical particles react with each other (for 
example in the 12 C + 12 C or the triple-ct reactions; for details see [59]). 

In addition to nuclear reactions, expansion or contraction of the plasma can 
also produce changes in the number densities rij. To separate the nuclear 
changes in composition from these hydrodynamic effects, we introduce the 
nuclear abundance Y^ = Hi/pN^, where Na is Avagadro's number. For a nu- 
cleus with atomic weight Aj, A^Yi represents the mass fraction of this nucleus, 
therefore J2AiYi = 1. Likewise, the equation of charge conservation becomes 
Y, ZiYi = Y e , where Y e (= n e /pN A ) is the electron abundance. By recasting 
Eq. 11 in terms of nuclear abundances Yi, a set of ordinary differential equa- 
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tions for the evolution of Yi results which depends only on nuclear reactions. In 
terms of the reaction cross sections introduced in §3.1, this reaction network 
is described by the following set of differential equations 

Y t = £ A/^Y} + J2KkP N A(J, QYjYk 

+ T,^kiP 2 N 2 A (j,k,l)Y ] Y k Y l . (12) 

j,k,l 



3. 3 Coupling Nuclear Networks to Hydrodynamics 



As we touched on in the previous section, nuclear processes are tightly linked 
to the hydrodynamic behavior of the bulk medium. Thermonuclear processes 
release (or absorb) energy, altering the pressure and causing hydrodynamic 
motions. These motions may disperse the thermonuclear ash, bringing a con- 
tinued supply of fuel to support the flame. The compositional changes, both 
of nuclei and of leptons, caused by thermonuclear reactions can also change 
the equation of state and opacity, further impacting the hydrodynamic behav- 
ior. For purposes of this review of nucleosynthesis methods, which generally 
assume that thermonuclear and hydrodynamic changes in local composition 
can be successfully decoupled, we include a brief description of how this de- 
coupling is best achieved. Miiller [109] provides an authoritative overview and 
discusses the difficulties (and open issues) involved when including nucleosyn- 
thesis within hydrodynamic simulations. 

The coupling between thermonuclear processes and hydrodynamic changes can 
be divided into two categories by considering the spatial extent of the cou- 
pling. Nucleosynthetic changes in composition and the resultant energy release 
produce local changes in hydrodynamic quantities like pressure and temper- 
ature. The strongest of these local couplings is the release (or absorption) of 
energy and the resultant change in temperature. Changes in temperature are 
particularly important because of the exponential nature of the temperature 
dependence of thermonuclear reaction rates. Since the nuclear energy release 
is uniquely determined by the abundance changes, the rate of thermonuclear 
energy release, e, is given by 

e n uc= -J2 N AM t c 2 Y l (MeV g- 1 s- 1 ). (13) 



where MiC 2 is the rest mass energy of species % in MeV. Since all reactions 
conserve nucleon number, the atomic mass excess M ex ^ = Mi — Aiin u (m u 
is the atomic mass unit) can be used in place of the mass Mi in Eq. 13 (see 
[110] for a recent compilation of mass excesses). The use of atomic mass units 



13 



has the added benefit that electron conservation is correctly accounted for in 
the case of f3~ decays and e~ captures, though reactions involving positrons 
require special treatment. In general, the nuclear energy release is deposited 
locally, so the rate of thermonuclear energy release is equal to the nuclear 
portion of the hydrodynamic heating rate. However, there are instances where 
nuclear products do not deposit their energy locally. Escaping neutrinos can 
carry away a portion of the thermonuclear energy release. In the rarefied envi- 
ronment of supernova ejecta at late times, positrons and gamma rays released 
by j3 decays are not completely trapped. In most such cases, the escaping par- 
ticles stream freely from the reaction site, allowing adoption of a simple loss 
term analogous to Eq. 13 with MiC 2 replaced by an averaged energy loss term. 
For example, the weak reaction rate tabulations of Fuller, Fowler, & Neumann 
[85] provide averaged neutrino losses. From these we can construct 

€v loss / \ *- J v) -* i,weakj v / 



where we consider only those contributions to Y due to neutrino producing 
reactions. In some cases, like supernova core collapse (see [20], this volume), 
more complete transport of the escaping leptons or gamma rays must be con- 
sidered. Other important quantities which are impacted by nucleosynthesis, 
like Y e , can be obtained by appropriate sums over the abundances and also 
need not be evolved separately. 

Implicit solution methods require the calculation of Y(t + At), where At is 
the nuclear timestep, which in turn requires knowledge of T(t + At). One 
could write a differential equation for the energy release analogous to Eq. 12, 
with the Afs replaced by the reaction Q-values, and thereby evolve the energy 
release (and calculate temperature changes) as an additional equation within 
the network solution. Muller [111] has shown that such a scheme can help avoid 
instabilities in the case of a physically isolated zone entering or leaving nuclear 
statistical equilibrium. In general, however, use of this additional equation 
is made unnecessary by the relative slowness with which the temperature 
changes. The timescale on which the temperature changes is given by 

t t = T/f « C v T/e nuc (15) 



and is often called the ignition timescale. The timescale on which an individual 
abundance changes is its burning time, 

r( A Z) = F( A Z)/F( A Z) = minr fe ( A Z) (16) 

k 



where Tfc( A Z) is defined in Eq. 4. In general tt differs from t( a Z) of the 
principle fuel by the ratio of thermal energy content to the energy released by 
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the reaction. For degenerate matter this ratio can approach zero, allowing for 
explosive burning. In contrast, best results for the nuclear network are achieved 
[112] when the network timestep At is chosen to be the burning timescale of 
a less abundant species, typically with an abundance of 10 -6 or smaller. Since 
the dominant fuel is typically one of the more abundant constituents and the 
burning timescales are proportional to the abundance, tt is typically an order 
of magnitude or more larger than the network timestep (see, e.g., [113,114]. It 
is therefore sufficient to calculate the energy gain at the end of a timestep via 
equation 13, modified as discussed above, and approximate T(t + St) ~ T(t) 
or to extrapolate based on e(t). Since other locally coupled quantities have 
characteristic timescales much longer than At, they too can be decoupled 
in a similar fashion. For the remainder of this review, we will consider only 
the equations governing changes in isotopic abundances, remembering that 
additional equations can easily be constructed for those special circumstances 
where they are necessary. 

Spatial coupling, particularly the modification of the composition by hydrody- 
namic movements such as diffusion, convective mixing and advection (in the 
case of Eulerian hydrodynamics methods), represents a more difficult chal- 
lenge. By necessity, an individual nucleosynthesis calculation examines the 
abundance changes in a locality of uniform composition. The difficulties as- 
sociated with strong spatial coupling of the composition occur because this 
nucleosynthetic calculation is spread over an entire hydrodynamic zone. Con- 
vection can result in strong abundance gradients across a single hydrodynamic 
zone, which with the assumption of compositional uniformity, can result in 
very different outcomes as a function of the fineness of the hydrodynamic 
grid. Eulerian advection of compositional boundaries can also have extremely 
unphysical consequences. Fryxell et al. [115] demonstrated how this artificial 
mixing can produce an unphysical detonation in a shock tube calculation by 
mixing cold unburnt fuel into the hot burnt region. A related problem is the 
conservation of species. Hydrodynamic schemes must carefully conserve the 
abundances (or partial densities) of all species [115-117], lest they provide un- 
physical abundances to the nucleosynthesis calculations, which must assume 
conservation, and thereby produce unphysical results. Because of these prob- 
lems, nucleosynthesis calculations are best suited to hydrodynamic simulations 
with excellent capture of shock and contact discontinuities. 

The relative size of the burning timescales, when compared to the relevant 
diffusion, sound crossing or convective timescale, dictates how these problems 
must be addressed. If all of the burning timescales are much shorter than 
the timescale on which the hydrodynamics changes the composition, then 
the assumption of uniform composition is satisfied and the nucleosynthesis 
of each hydrodynamic zone can be treated independently. If all of the burn- 
ing timescales are much larger than, for example, the convective timescale, 
then the composition of the entire convective zone can be treated as uni- 
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form and slowly evolving. The greatest complexity occurs when the timescales 
on which the hydrodynamics and nucleosynthesis change the composition are 
similar. Oxygen shell burning represents an excellent example of this as the 
sound travel, convective turnover and nuclear burning timescales are all of 
the same order as the evolutionary time. The results (of 2D simulations [118]) 
demonstrate convective overshooting, highly non-uniform burning and a veloc- 
ity structure dominated by convective plumes. Silicon burning also represents 
a particular challenge [7], as the timescales for the transformation of silicon 
to iron are much slower than the convective turnover time, but the burning 
timescales for the free neutrons, protons and a-particles which maintain QSE 
are much faster, providing a strong motivation for the hybrid networks we will 
discuss in 86. 



4 Solving the Nuclear Network 

In principle, the initial value problem presented by the nuclear network can be 
solved by any of a large number of methods discussed in the literature. However 
the physical nature of the problem, reflected in the A's and (cru)'s, greatly 
restricts the optimal choice. The large number of reactions display a wide range 
of reaction timescales, r (see Eq. 4). Systems whose solutions depend on a wide 
range of timescales are termed stiff. Gear [119] demonstrated that even a single 
equation can be stiff if it has both rapidly and slowly varying components. 
Practically, stiffness occurs when the limitation of the timestep size is due to 
numerical stability rather than accuracy. A more rigorous definition [120] is 

that a system of equations Y(Y) is stiff if the eigenvalues Xj of the Jacobian 

dY/dY obey the criteria 

»(A J -)<0, j = l,---,N (17) 

s _ max\^( y \ j )\ i 
min 1 9ft(Aj) | 

where Sft(A) is the real part of the eigenvalues A. As we will explain in this 
section, S > 10 15 is not uncommon in astrophysics. 

Nucleosynthesis calculations belong to the more general field of reactive flows, 
and therefore share some characteristics with related terrestrial fields. In par- 
ticular, chemical kinetics, the study of the evolution of chemical abundances, 
is an important part of atmospheric and combustion physics and produces 
sets of equations much like Eq. 11 (see [121] for a good introduction). These 
chemical kinetics systems are known for their stiffness and a great deal of ef- 
fort has been expended on developing methods to solve these equations. Many 
of the considerations for the choice of solution method for chemical kinetics 
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also apply to nucleosynthesis calculations. In both cases, temporal integration 
of the reaction rate equations is broken up into short intervals because of the 
need to update the hydrodynamics variables. This favors one step, self starting 
algorithms. Because abundances must be tracked for a large number of com- 
putational cells (hundreds to thousands for one dimensional models, millions 
for the coming generation of three dimensional models), memory storage con- 
cerns favor low order methods since they don't require the storage of as much 
data from prior steps. In any event, both the errors in fluid dynamics and in 
the reaction rates are typically a few percent or more, so the greater precision 
of these higher order methods often does not result in greater accuracy. 

Because of the wide range in timescales between strong, electromagnetic and 
weak reactions, the nuclear networks are extraordinarily stiff. PP chain nu- 
cleosynthesis, responsible for the energy output of the Sun, offers an excellent 
example of the difficulties. The first reaction of the PP1 chain is 1 H(p, e + z/) 2 H, 
the fusion of two protons to form deuterium. This is a weak reaction, requiring 
the conversion of a proton into a neutron, and releasing a positron and a neu- 
trino. As a result, the reaction timescale r p ( 1 H) is very long, billions of years 
for conditions like those in the solar interior. The second reaction of the PP1 
chain is the capture of a proton on the newly formed deuteron, 2 H(p, 7) 3 He. 
For conditions like those in the solar interior, the characteristic timescale, 
r p ( 2 H) is a few seconds. Thus the timescales for two of the most important 
reactions for hydrogen burning in stars like our Sun differ by more than 17 
orders of magnitude (see [5] for a more complete discussion of the PP chain). 
This disparity results not from a lack of p+p collisions (which occur at a rate 
y( 1 H)/y( 2 H) ~ 10 17 times more often than 1 H + 2 H collisions), but from 
the rarity of the (weak) transformation of a proton to a neutron. While the 
presence of weak reactions among the dominant energy producing reactions is 
unique to hydrogen burning, most nucleosynthesis calculations are similarly 
stiff, in part because of the need to include weak interactions but also the 
potential for neutron capture reactions, which occur very rapidly even at low 
temperature, following any release of free neutrons. Though further investi- 
gation is warranted, the nature of the nuclear reaction network equations has 
thus far limited the astrophysical usefulness of the most sophisticated methods 
to solve stiff equations developed for chemical kinetics. 

For a set of nuclear abundances Y, one can calculate the time derivatives of 

the abundances, Y using Eq. 12. The desired solution is the abundance at a 
future time, Y(t + At), where At is the network timestep. Since coupling with 
hydrodynamics favors low order, one step methods, general nucleosynthesis 
calculations use the simple finite difference prescription 

Y{t + A £~ Y{t) = (1 - 0)f(t + At) + 0f(t). (18) 
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With = 1, Eq. 18 becomes the explicit Euler method while for O = it is the 
implicit backward Euler method, both of which are first order accurate. For 
= 1/2, Eq. 18 is the semi-implicit trapezoidal method, which is second order 
accurate. For the stiff set of non-linear differential equations which form most 
nuclear networks, a fully implicit treatment is generally most successful [112], 
though the semi-implicit method has been used in Big Bang nucleosynthesis 
calculations [52], where coupling to hydrodynamics is less important. Solving 
the fully implicit version of Eq. 18 is equivalent to finding the zeros of the set 
of equations 

This is done using the Newton- Raphson method (see, e.g., [122]), which is 
based on the Taylor series expansion of Z(t + At), with the trial change in 
abundances given by 

AY = (»$t±M)" g , (20) 



where dZ/dY is the Jacobian of Z. Iteration continues until Y(t + At) con- 
verges. 

A potential numerical problem with the solution of Eq. 19 is the singularity 
of the Jacobian matrix, dZ (t + At)/dY(t + At). From Eq. 19, the individual 
matrix elements of the Jacobian have the form 



dZ t 
dY 3 


At 


dY % 
dY 3 




At 


-E l 



(21) 



where 5^ is the Kronecker delta, and Tj{i) is the destruction timescale of nu- 
cleus % with respect to nucleus j for a given reaction, as defined in Eq. 4. The 
sum accounts for the fact that there may be more than one reaction by which 
nucleus j is involved in the creation or destruction of nucleus %. Along the 
diagonal of the Jacobian, there are two competing terms, 1/At and J2 1/ t j(0- 
This sum is over all reactions which destroy nucleus i, and is dominated by the 
fastest reactions. As a result, J2 l/ r i W can be orders of magnitude larger than 
the reciprocal of the desired timestep, 1/At. This is especially a problem near 
equilibrium, where both destruction and the balancing production timescales 
are very short in comparison to the preferred timestep size, resulting in differ- 
ences close to the numerical accuracy (i.e. 14 or more orders of magnitude). In 



such cases, the term 1/At is numerically neglected, leading to numerically sin- 
gular matrices. One approach to avoiding this problem is to artificially scale 
these short, equilibrium timescales by a factor which brings their timescale 
closer to At, but leaves them small enough to ensure equilibrium. While this 
approach has been used successfully, the ad hoc nature of this artificial scaling 
renders these methods fragile. A more promising approach is to make directly 
use of equilibrium expressions for abundances, which, as we will discuss in §6, 
also assures the economical use of computer resources. 



4-1 Taking Advantage of Matrix Sparseness 



For larger networks, the Newton- Raphson method requires solution of a mod- 
erately large (N = 100 — 3000) matrix equation. Since general solution of a 
dense matrix scales as 0(N 3 ), this can make these large networks progres- 
sively much more expensive. While in principal, every species reacts with each 
of the hundreds of others, resulting in a dense Jacobian matrix, in practice it is 
possible to neglect most of these reactions. Because of the Z t Zj dependence of 
the repulsive Coulomb term in the nuclear potential, captures of free neutrons 
and isotopes of H and He on heavy nuclei occur much faster than fusions of 
heavier nuclei. Furthermore, with the exception of the Big Bang nucleosyn- 
thesis and PP-chains, reactions involving secondary isotopes of H (deuterium 
and tritium) and He are neglectable. Likewise, photodisintegrations tend to 
eject free nucleons or a-particles. Thus, with a few important exceptions, for 
each nucleus we need only consider twelve reactions linking it to its nuclear 
neighbors by the capture of an n, p, a or 7 and release a different one of these 
four. The exceptions to this rule are the few heavy ion reactions important 
for burning stages like carbon and oxygen burning where the dearth of light 
nuclei cause the heavy ion collisions to dominate. 

Fig. 2 demonstrates the sparseness of the resulting Jacobian matrix, for a 300 
nuclei network designed for silicon burning, but capable of handling all prior 
burning stages. Of the 90,000 matrix elements, less than 5,000 are non-zero. In 
terms of the standard forms for sparse matrices, this Jacobian is best described 
as doubly bordered, band diagonal. With a border width, AB, of 45 necessary 
to include the heavy ion reactions among 12 C, 16 and 20 Ne along with the 
free neutrons, protons and a-particles and a band diagonal width, AD, of 
54, even this sparse form includes almost 50,000 elements. With solution of 
the matrix equation consuming 90+% of the computational time, there is 
clearly a need for custom tailored solvers which take better advantage of the 
sparseness of the Jacobian. To date best results for small (N<100) matrices 
are obtained with machine optimized dense solvers (e.g. LAPACK) or matrix 
specific solvers generated by symbolic processing [111,109]. For large matrices, 
generalized sparse solvers, both custom built and from software libraries, are 
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Fig. 2. Graphic demonstration of the sparseness of the Jacobian matrix. The filled 
squares represent the non-zero elements. 

used (see, e.g., [123]). 



4-2 Physically Motivated Network Specialization 



Often from a physical understanding one can specialize the general solution 
method and thereby greatly reduce the computational cost. As an example of 
such, we will in this section discuss the r-process approximation of [73] (see 
also [124,125]). For nuclei with A > 100, charged particle captures (proton 
and a) as well as their reverse photodisintegrations virtually cease when T < 
3 GK. This leaves only neutron captures and their reverse photodisintegration 
reactions, as well as /3-decays, which can also lead to the emission of delayed 
neutrons (here we consider the release of up to three delayed neutrons). In 
this case, Eq. 12 greatly simplifies, leaving 
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y( A z) = n n (^)^_ 1 y( A - 1 z) + a^ +1 f( a+1 z) + £ Ag_V+^( A+j z - 1) 

nnMSJ + A1, A + E A?3 ) y ( A Z) , (22) 



where (o"f)^'3i and A^ are the velocity integrated neutron capture cross sec- 
tion and the photodisintegration rate for the nucleus A Z, while X^^A) * s ^ ne 
decay constant for the f3~ decay of A Z, with j delayed neutrons. The assump- 
tion is made that the neutron abundance (Y n = n n / pNjC) varies slowly enough 
that it may be evolved explicitly. One can see that in Eq. 22, with n n thereby 
fixed, the time derivatives of each species have a linear dependence on only 
the abundances of their neighbors in the same isotopic chain (nuclei with the 
same Z), or that with one less proton (Z — 1). One can then divide the network 
into separate pieces for each isotopic chain, and solve them sequentially, be- 
ginning with the lowest Z. The "boundary" terms for this lowest Z chain can 
be supplied by a previously run or concurrently running full network calcula- 
tion which need extend only to this Z. This reduces the solution of a matrix 
with more than a thousand rows to the solution of roughly 30 smaller matrices. 
Furthermore each of these smaller matrices is also tridiagonal increasing speed 
further. Freiburghaus et al. [125] tested the assumption of slow variation in 
the neutron abundance, and have demonstrated the usefulness of this method 
in r-process simulations, achieving a large decrease in computational cost. A 
similar treatment has been successfully applied to explosive hydrogen burn- 
ing based on the assumption of slowly varying proton and alpha abundances 
[47]. As we will discuss in §6, for other burning stages there exist physically 
motivated simplifications to the general network solution method. 



5 Equilibria in Nuclear Astrophysics 



As is the case in many disciplines, equilibrium expressions are frequently em- 
ployed to simplify nuclear abundance calculations. In most such cases of inter- 
est in nuclear astrophysics, the fast strong and electromagnetic reactions reach 
equilibrium while those involving the weak nuclear force do not. Since the weak 
reactions are not equilibrated, the resulting Nuclear Statistical Equilibrium 
(NSE) requires monitoring of weak reaction activity. Even with this stric- 
ture, NSE offers many advantages, since hundreds of abundances are uniquely 
defined by the thermodynamic conditions and a single measure of the weak 
interaction history or the degree of neutronization. Computationally, this re- 
duction in the number of independent variables greatly reduces the cost of nu- 
clear abundance evolution. Because there are fewer variables to follow within 
a hydrodynamic model, the memory footprint of the nuclear abundances is 
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also reduced, an issue of importance in modern multi-dimensional models of 
supernovae. Finally, the equilibrium abundance calculations depend on bind- 
ing energies and partition functions, quantities which are better known than 
many reaction rates. This is particularly true for unstable nuclei and for con- 
ditions where the mass density approaches that of the nucleus itself, resulting 
in exotic nuclear structures. 

The expression for NSE is commonly derived using either chemical potentials 
or detailed balance (see, e.g., [5,126-129]). For a nucleus A Z, composed of Z 
protons and N = (A—Z) neutrons, in equilibrium with these free nucleons, the 
chemical potential of A Z can be expressed in terms of the chemical potentials 
of the free nucleons 

Hz,a = Zfi p + Nfx n . (23) 



For a collection of particles obeying Boltzmann statistics, the chemical poten- 
tial, including rest mass, of each species is given by 
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(24) 



(e.g., [130]). Substituting Eq. 24 into Eq. 23 allows derivation of an expression 
for the abundance of every nuclear species in terms of the abundances of the 
free protons (Y p ) and neutrons (Y n ), 



,A*s G ( AZ )(P N A\ A -\l„(B{ A Z)\ v „ 



= C( A Z)Y n N Y p z , (25) 

where G( A Z) and B( A Z) are the partition function and binding energy of the 
nucleus A Z, Na is Avagadro's number, k B is Boltzmann's constant, p and T 
are the density and temperature of the plasma, and 9 is given by 



' m u k B T \ 
2nh 2 I 



(26) 



Thus abundances of all nuclear species can be expressed as functions of two. 
Mass conservation (J2 AY = 1) provides one constraint. The second constraint 
is the amount of weak reaction activity, often expressed in terms of the total 
proton abundance, J2 ZY, which charge conservation requires equal the elec- 
tron abundance, Y e . Thus the nuclear abundances are uniquely determined for 
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Fig. 3. The average atomic mass for material in NSE as a function of Temperature. 
The solid lines include screening corrections to the nuclear binding energies, while 
the dotted lines ignore this effect. See [129] for more discussion of the importance 
of screening in NSE. 

a given (T,p,Y e ). Alternately, the weak interaction history is sometimes ex- 
pressed in terms of the neutron excess 77 = XX iV — Z)Y . Figure 3 displays the 
temperature and density dependence of A = X AY/ Y^Y — 1/Y^Y, the aver- 
age nuclear mass of the NSE distribution. At high temperatures, free nucleons 
are favored, hence A ~ 1. For intermediate temperatures the compromise of 
retaining large numbers of particles while increasing binding energy favors 4 He, 
which has 80% of the binding energy of the iron peak nuclei. At low temper- 
atures, Eq. 25 strongly favors the most bound nuclei, the iron peak nuclei, so 
A — > 60 as the temperature drops. Density can be seen to scale the placement 
of these divisions between high, intermediate and low temperature. Variations 
in Y e do not strongly affect Figure 3. At high temperatures, it simply effects 
the ratio of Y p /Y n rs Y e /1 — Y e . At low temperatures, variation in Y e changes 
which Fe-peak isotopes dominate. Though 56 Ni is less tightly bound than 54 Fe, 
it is more tightly bound than 54 Fe + 2 1 H, which would be required by charge 
conservation if Y e ~ .5. Thus F( 56 Ni) > F( 54 Fe) for low T with Y e ~ .50, but 
y( 54 Fe) > y( 56 Ni) for smaller Y e . In general, the most abundant nuclei at low 
temperatures are the most bound nuclei for which Z/A ~ Y e . 

As with any equilibrium distribution, there are limitations on the applica- 
bility of NSE. The first requirement for NSE to provide a good estimate of 
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the nuclear abundances is that the temperature be sufficient for the endoer- 
gic reaction of each reaction pair to occur. Since for all particle-stable nuclei 
between the proton and neutron drip lines (with the exception of nuclei unsta- 
ble against alpha decay), the photodisintegrations are endoergic, with typical 
Q- values among ((3) stable nuclei of 8-12 MeV, by Eq. 9 this requirement 
reduces to T > 3GK. While this requirement is necessary, it is not sufficient. 
In the case of hydrostatic silicon burning, even when this condition is met, 
appreciable time is required to convert Si to Fe-peak elements. In the case of 
explosive silicon burning, the adiabatic cooling on timescales of seconds can 
cause conditions to change more rapidly than NSE can follow, breaking NSE 
down first between 4 He and 12 C, at T ~ 6GK [131] and later between the 
species near silicon and the Fe-peak nuclei, at T ~ 4GK [132]. Thus it is clear 
that in the face of sufficiently rapid thermodynamic variations, NSE provides 
a problematic estimate of abundances. 



6 Merging Equilibria with Nuclear Networks 



In spite of the limitations on the applicability of NSE, the reduced compu- 
tational cost provides a strong motivation to maximize the use of equilibria. 
The use of equilibrium expressions for single abundances is, in fact, common in 
nuclear reaction networks, typically to track the abundances of short-lived un- 
stable intermediates in "three-particle" processes. The most common example 
of this is the triple a process, 



4 He + 4 He ^ 8 Be Q = -.09 MeV 

(27) 

4 He + 8 Be # 12 C* -> 12 C + 7 Q = 7.37 MeV , 



by which Helium burning occurs. With r( 8 Be) ~ 10~ 16 sec, only rarely does 
a 8 Be survive long enough for a second a to capture. As a result of the near 
balance of the first reaction pair, the abundance of 8 Be can be expressed in 
terms of the a-particle abundance, 



Likewise for temperatures in excess of . 1 GK, the most likely result following 
the second a capture to form an excited state of 12 C is a decay back to 8 Be 
(r a ( 12 C*)/r 7 ( 12 C*) > 10 3 ), thus the abundance of 12 C* is well characterized, 
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via/i( 12 C*) = 3p( 4 He),by 
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When this is the case, the effective triple a reaction rate is simply that of the 
decay of 12 C from the excited state to the ground state, 

r 3a = P iv A r( 12 c*)r 7 ( 12 c*)// i . (30) 



This use of local equilibrium within a rate equation shares many characteristics 
with the more elaborate schemes we will discuss later in this section. The 
number of species tracked by the network is reduced since Y( 8 Be) need not 
be directly evolved. Problematically small timescales like r( 8 Be) are removed, 
replaced by larger time scales {r 3a ~ 10 5 — 10 7 years during core helium 
burning). The non-linearity of network time derivatives is increased (Y( 12 C) oc 
Y£) under this scheme. This approximation also breaks down at low T (for 
details see [107]). 

In addition to silicon burning, there are a number of astrophysically impor- 
tant situations where T > Q/30kB for at least some of the relevant reactions 
and so large equilibrium groups exist, but NSE is not globally valid. This in- 
clude the r-process [73,133,134] and the rp-process in novae and X-ray bursts 
[47,51], where neutron or proton separation energies (Q n or Q p ) of 2 MeVand 
less are often encountered. Beginning with [135], a number of attempts have 
been made to take advantage of these partial equilibria to reduce the num- 
ber of independent variables evolved via rate equations and thereby reduce 
the computational cost of modeling these burning stages. Since we lack the 
time to discuss all in detail, we refer the reader to [36] and [47] for discussion 
of hybrid networks for the r-process and rp-process, respectively. Instead, we 
will here concentrate on the application of hybrid equilibria networks to sil- 
icon burning [11,131,132,136]. In the present context we will concentrate on 
the QSE-reduced a-network [137], a simple, but pedagogically illustrative, ex- 
ample which details the main ideas behind such hybrid equilibrium networks. 
Quasi-equilibrium (QSE) is a term coined by Bodansky, Clayton & Fowler 
[135] to describe the local equilibrium groups which form during silicon burn- 
ing. 



6.1 The QSE-reduced a-network 



Tracking the nuclear evolution during the major energy producing burning 
stages from the exhaustion of hydrogen through to the establishment of NSE 
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requires, at minimum, a network that includes nuclei from a-particles to Zn. 
As we discussed in §2, silicon burning presents a particular problem as material 
proceeds from silicon to the iron peak not via heavy ion captures but through a 
chain of photodisintegrations and light particle captures. We will discuss here 
the minimal nuclear set which can follow this evolution, the set of a-particle 
nuclei; a, 12 C, 16 0, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 52 Fe, 56 Ni, 
60 Zn. For convenience we will label this full set JF and refer to its abundances 
as Y T . Silicon burning in fact presents a larger problem, as the nuclear flow 
from silicon to the iron peak nuclei does not generally proceed through nuclei 
with N=Z, especially when significant neutronization has occurred [11]. In 
some hydrodynamical models, however, such compromise is made necessary 
by the computational limitations, either the time necessary to solve larger 
networks or the hydrodynamical problems associated with evolving and storing 
a large number of abundances. Furthermore, the small size of the a network (14 
nuclei and 17 reactions) makes application of QSE to a-chain nucleosynthesis 
a pedagogically useful example. 

The objective of the QSE-reduced a-network is to evolve Y T (and calculate 
the resulting energy generation) in a more efficient way. Under conditions 
where QSE applies, the existence of the silicon and iron peak QSE groups 
(which are separated by the nuclear shell closures Z=N=20 and the resulting 
small Q-values and reaction rates) allows calculation of these 14 abundances 
from 7. For the members of the silicon group ( 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti) and 
the iron peak group ( 48 Cr, 52 Fe, 56 Ni, 60 Zn) the individual abundances can be 
calculated by expressions similar to Eq. 25, 



C( A Z) A=2§. 

C( A Z) A-56 

Y QS eM A Z) = 7 ^^y( 56 Ni)F— , (31) 



C( 56 Ni) 

where C( Z) is defined in Eq. 25 and {A— 28)/4 and {A— 56)/4 are the number 
of a-particles needed to construct A Z from 28 Si and 56 Ni, respectively. Thus, 
where QSE applies, Y T is a function of Y n , where the reduced nuclear set 1Z 
is defined as a, 12 C, 16 0, 20 Ne, 24 Mg, 28 Si, 56 Ni, and we need only evolve Y n . 
It should be noted that 24 Mg is ordinarily a member of the silicon QSE group 
[7,10,11], but for easier integration of prior burning stages with a conventional 
nuclear network, we will evolve 24 Mg independently. The main task when 
applying such hybrid schemes is finding the boundaries of QSE groups and 
where individual nuclei have to be used instead. Treating marginal group 
members as part of a group increases the efficiency of the calculation, but 
may decease the accuracy. 

While Y 11 is a convenient set of abundances for calculating Y^, it is not the 
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most efficient set to evolve, primarily because of the non-linear dependence 
on Y a . Instead we define Y g = [Y aG , F( 12 C), F( 16 0), F( 20 Ne), F( 24 Mg), Y SiG , 
Y FeG ] where 



i£Si group iaFe group 

Y SiG = E Y, ( 32 ) 

idSi group 

^FeG — E Y% ■ 

iGFe group 

Physically, Y aG represents the sum of the abundances of free a-particles and 
those a-particles required to build the members of the QSE groups from 28 Si 
or 56 Ni, while Y$i G and Yp eG represent the total abundances of the silicon 
and iron peak QSE groups. This method, which here is applied only to the 
chain of a-nuclei can also be generalized to arbitrary networks [136]. For larger 
networks which contain nuclei with N ^ Z, one must be able to follow the 
abundances of free neutrons and protons, particularly since weak interactions 
will change the global ratio of neutrons to protons. In place of Y aG in Eq. 32, 
one constructs Y NG = 52i,u g ht N i Y i + T,i,Si( N i ~ l4 ) Y i + Y,i,Fe( N i ~ 28)3^ and 
Y ZG = Y.i,u gh t Z t Y t + Ei, Si ( Z i - U )Y + Ei, Fe ( Z i ~ 2S )Y, if 28 Si and 56 Ni are 
chosen as the focal nuclei for the Si and Fe groups. 

Corresponding to this reduced set of abundances Q is a reduced set of reac- 
tions, with quasi-equilibrium allowing one to ignore the reactions among the 
members of the QSE groups. Unfortunately, the rates of these remaining reac- 
tions are functions of the full abundance set, Y^, and are not easily expressed 
in terms of the group abundances, Y* 5 . Thus, for each Y G , one must solve for 

Y n and, by Eq. 31, Y^, in order to calculate Y* 3 which is needed to evolve Y* 5 
via Eq. 18. Furthermore, Eq. 20 requires the calculation of the Jacobian of Z, 

which can not be calculated directly since Y$ cannot be expressed in terms of 
Y g . Instead we find it sufficient to use the chain rule, 

QyQ _ QfQ dY n 
W^ ~ dY n W~s 



to calculate the Jacobian. Analytically, the first term of the chain rule product 
is easily calculated from the sums of reaction terms, while the second term 
requires implicit differentiation using Eq. 32, which is discussed further by Hix 
et al. [137]. 
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6.2 Silicon burning with the QSE-reduced a network 



In this section we will demonstrate the accuracy with which the QSE-reduced 
a network duplicates the results of the full 14 element a network for silicon 
burning. Our first example is a nucleosynthesis calculation occurring under 
constant temperature and density. While such calculations provide the least 
challenging comparison, they also allow comparison with NSE, which should 
represent the final abundances for these calculations. Figure 4 offers com- 
parison of the mass fractions of the 7 independent species; a-particles, 12 C, 
16 0, 20 Ne, 24 Mg, and the silicon and iron peak groups, as evolved by the QSE- 
reduced and conventional a networks for silicon burning at 5 GK and a density 
of 10 9 gcm~ 3 . Apart from an early enhancement by the QSE-reduced network 
of the iron peak mass fraction (20% after 10~ 6 seconds), these mass fractions 
typically agree to within 1%. Since the nuclear energy release depends linearly 
on the abundance changes (see Eq. 13), differences in small abundances have 
little effect on the nuclear energy store. In this case, the difference in the rate 
of energy generation calculated by the two networks is < 1% at 10~ 6 and 
10~ 4 seconds. This difference is significantly smaller than the variation, shown 
by both networks, in the rate of energy generation between timesteps, with e 
typically declining by 5% per timestep over this interval. 

Significant variations in abundance among the individual members of the QSE 
groups between the two networks are also limited to the small abundances. At 
early times, the small abundances within the iron peak reduce the accuracy 
of QSE at predicting the individual abundances of members of the iron peak 
group. Much of the enhanced mass fraction of the iron peak nuclei at early 
times, seen in Fig. 4, is due to the QSE-reduced network's emphasis on heavier 
nuclei at the expense of 48 Cr. After an elapsed time of 10~ 6 seconds, the average 
mass of the iron peak nuclei, Ap e G, is 49.2 according to the conventional 
network and 52.6 according to the QSE-reduced network. As a result, the 
abundances of 48 Cr and 52 Fe calculated by the QSE-reduced network are 38% 
and 164% of their conventional network values, while 56 Ni and 60 Zn are 16 
times more abundant than the conventional network predicts. However, the 
total Fe group mass fraction is only 10 -7 at this point in time. As the iron 
peak nuclei become more abundant, QSE provides a better estimate of the 
relative abundances within the group, reducing such discrepancies. By the 
time the iron peak nuclei represent a significant portion of the mass, the 
differences in the abundance predictions for all nuclei are only a few per cent. 
As each network reaches its respective equilibrium, after an elapsed time of 
24 seconds, the abundance predictions of these networks (shown in Table 2) 
differ by at most 3%, even among the nuclei with the smallest abundances. Not 
surprisingly, in view of these small abundance differences, the difference in the 
total energy released by these networks is less than .1%. Comparison of the 
network abundances with abundances calculated from NSE show a similarly 
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Fig. 4. Evolution of the independent nuclear mass fractions for constant thermo- 
dynamic conditions, T = 5 GK and p = 10 9 gcm~ 3 . The solid lines display the 
evolution due to a conventional a-network, the circles show the evolution by the 
QSE-reduced a network. 

low level of difference. 

While example calculations of silicon burning under constant conditions are in- 
structive, if the QSE-reduced a-network is to replace a conventional a-network 
it must be shown to be accurate under changing thermodynamic conditions. 
Of particular importance is the ability to model explosive silicon burning. 
Because the products of hydrostatic silicon burning are trapped deep in the 
potential well of their parent star, it is only by explosion that the interstellar 
medium is enriched in intermediate mass and iron peak elements. To model 
silicon burning occurring as a result of shock heating, we will follow the approx- 
imation introduced by [138]. Therein a mass zone is instantaneously heated by 
a passing shock to some peak initial temperature, T t , and density, p iy and then 
expands and cools adiabatically, with the timescale for the expansion given 



by the free fall timescale, thd = (24irGp) 



-1/2 



— 1/2 

446p 6 milliseconds. 



Figure 5 shows the nuclear evolution for an example of this explosive burning 
model with Tj = 5GK and pi = 10 9 gcm~ 3 . Over the first millisecond, the 
evolution portrayed here is virtually identical to that of Figure 4; however by 
the time one millisecond has elapsed, the temperature has dropped to 4.9 GK, 
slowing the reactions which are turning silicon into iron peak elements. This 
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Table 2 

Comparison of equilibria calculated by the conventional and QSE-reduced a- 

networks with NSE for T = 5 GK and p = 10 9 g cm~ 3 



Nucleus 


*net 


Y 

1 qse 


Y 

1 nse 


4 He 


7.73 xl(T 5 


7.78xl(T 5 


7.80xl(T 5 


12 C 


l.OlxlO" 10 


1.02xl0~ 10 


l.OlxlO -10 


16Q 


3.23xl(T 10 


3.30xl(T 10 


3.24xl(T 10 


20 Ne 


4.08xl(T 12 


4.19xl(T 12 


4.11xl(T 12 


24 Mg 


2.13xl(T 9 


2.20xl(T 9 


2.16xl(T 9 


28 Si 


4.51xl(T 6 


4.68xl(T 6 


4.53xl(T 6 


32 S 


9.34xl(T 6 


9.63xl(T 6 


9.37xl(T 6 


36 Ar 


1.06xl(T 5 


1.09xl(T 5 


1.07xl(T 5 


40 Ca 


3.01xl(T 5 


3.06xl(T 5 


3.03xl(T 5 


44 Ti 


1.68xl(T 6 


1.70xl(T 6 


1.69xl0 -6 


48 Cr 


3.75xl(T 5 


3.80xl(T 5 


3.75xl(T 5 


52 Fe 


8.64xl(T 4 


8.70xl(T 4 


8.65xl(T 4 


56 M 


1.70xl0 -2 


1.70xl(T 2 


1.70xl(T 2 


60 Zn 


4.71xl(T 6 


4.68xl(T 6 


4.69xl(T 6 



cooling, which drops T below 4 GK after 9 ms and below 3 GK after 22 ms, 
freezes out the nuclear reactions before NSE is reached, resulting in incomplete 
silicon burning, as discussed by Woosley et al. [10]. In this case, the freezeout 
leaves nearly equal amounts of silicon group and iron peak group elements. 
Throughout most of the evolution in this example, the agreement between the 
mass fractions as evolved by the QSE-reduced a network with those evolved 
by its conventional counterpart is comparable to that demonstrated under 
constant thermodynamic conditions. Columns 2 and 3 of Table 3 compare the 
abundances after 9 ms have elapsed, with T nearing 4 GK. In this case, the 
largest relative error (5%) is in the abundance of 48 Cr. These small differences 
in abundance result in small differences in the accumulated nuclear energy 
release, approximately .5% to this point. By this time, adiabatic cooling has 
greatly reduced the rate of energy generation from its peak of more than 
10 22 ergsg _1 s^ 1 to roughly 10 17 ergsg _1 s" 1 . Though the absolute difference in 
the rate of energy generation as calculated by the two networks has declined 
from ~ 10 19 to 10 16 ergsg _1 s _1 , the relative difference has grown to 10% as 
T nears 4 GK. Fortunately this difference is negligible, since nuclear energy 
release has virtually ceased, with < .2% of the total energy release remaining. 

Though energetically unimportant, nuclear reactions continue, resulting in sig- 
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Fig. 5. Evolution of the independent nuclear mass fractions under adiabatic expan- 
sion with Tj = 5 GK and pi = 10 9 gem -3 . The solid lines display the evolution due 
to a conventional a network, the circles show the evolution by the QSE-reduced a 
network. 

nificant changes in the smaller abundances, a point which is discussed in detail 
by [10,132]. The continued cooling drops T below 3 GK, gradually freezing 
out the photodisintegrations responsible for QSE. However this same decline 
in temperature also reduces the rate of charged particle capture reactions, 
greatly reducing the amount of nucleosynthesis which occurs after T drops 
below ~ 3.5 GK. The group abundances of the silicon and iron peak groups 
(which account for 99.9% of the mass), as calculated by the QSE-reduced a 
network after 18 ms have elapsed (T = 3.3 GK), differ by less than 1% from 
those of the conventional a network at the same point in time. Comparison 
of Columns 4 and 5 of Table 3 reveals larger variations among individual 
abundances, most notably, significant under estimation by the QSE-reduced 
network of the 48 Cr and 52 Fe abundances, with a small, compensatory over 
estimate of the 56 Ni abundance. These variations, factors of 2 and 4 for 52 Fe 
and 48 Cr, respectively, and 3% for 56 Ni signal the breakdown of the iron peak 
QSE group. With the steep decline in temperature and density, the flux up- 
ward from 52 Fe in the conventional network is no longer sufficient to provide 
the reduction in abundance which QSE and the sharply declining abundance 
of free a-particles requires. 

As the temperature and density continue to decline, so too does the free a- 
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Table 3 

Comparison of network abundances for Tj = 5.0 GK and pi= 10 9 gcm~ 3 



Time (ms) 




8.77 




17.7 


255 




T(GK) 




4.07 




3.29 


0.01 




Nucleus 


inet 




Y 

1 qse 




-'net 


Y 

1 qse 


* net 




4 He 


7.90x10" 


-7 


7.82x10" 


-7 


1.04xl0~ 8 


1.01xl0~ 8 


1.94x10" 


-14 


12 C 


1.96x10" 


-7 


1.96x10" 


-7 


3.99xl0~ 8 


3.23xl0~ 8 


3.90x10" 


-8 


16q 


7.34x10" 


-7 


7.39x10" 


-7 


5.26xl0~ 7 


5.07xl0~ 7 


5.27x10" 


-7 


20 Ne 


2.63x10" 


-9 


2.63x10" 


-9 


1.69xl0- 10 


1.48X10" 10 


9.88x10" 


-11 


24 Mg 


2.24x10" 


-6 


2.26x10" 


-6 


2.88xl0~ 7 


2.98xl0~ 7 


2.80x10" 


-7 


28 Si 


7.65x10" 


-3 


7.76x10" 


-3 


7.58xl0~ 3 


7.86xl0~ 3 


7.58x10" 


-3 


32 S 


4.93x10" 


-3 


4.96x10" 


-3 


5.16xl0~ 3 


5.15xl0~ 3 


5.16x10" 


-3 


36 Ar 


1.43x10" 


-3 


1.42x10" 


-3 


1.27xl0~ 3 


1.21xl0~ 3 


1.27x10" 


-3 


40 Ca 


1.32x10" 


-3 


1.30x10" 


-3 


1.32xl0~ 3 


1.22xl0~ 3 


1.32x10" 


-3 


44 Ti 


7.07x10" 


-6 


6.90x10" 


-6 


1.96xl0~ 6 


1.72xl0~ 6 


1.69x10" 


-6 


48 Cr 


5.89x10" 


-5 


5.58x10" 


-5 


4.40xl0~ 5 


1.19xl0~ 5 


4.40x10" 


-5 


52 Fe 


7.17x10" 


-4 


6.93x10" 


-4 


6.33xl0~ 4 


3.05xl0~ 4 


6.33x10" 


-4 


56 Ni 


8.63x10" 


-3 


8.60x10" 


-3 


8.73xl0~ 3 


9.04xl0~ 3 


8.73x10" 


-3 


60 Zn 


6.18x10" 


-8 


6.06x10" 


-8 


3.26xl0~ 9 


3.23xl0~ 9 


4.38x10" 


-10 



particle abundance. Column 6 of Table 3 details the abundances after 255 ms 
have elapsed, with T having dropped to .01 GK, and all abundances having 
frozen out. Comparison of columns 4 and 6 reveals that the decline of the free 
a-particle abundance from 10~ 8 is the largest abundance variation beyond 
18 ms. Since the more abundant species have effectively frozen out by the 
time T approaches 3.5 GK, comparison of columns 5 and 6 reveal that the 
predictions of the QSE-reduced network, frozen near T = 3.5 GK, also provide 
good abundance estimates, in spite of the effects of differential freezeout. Thus, 
we have seen that the QSE-reduced a network provides an excellent description 
of energy generation and an accurate account of the abundances, down to 
(reaction) freeze-out temperatures of ~ 3.5 GK. These abundances at T ~ 
3.5 GK provide a very good approximation to the final results for the dominant 
nuclei. More detailed and accurate accounting of smaller abundances would 
require a switch back to the use of the full network below these freeze-out 
temperatures. 

In this section, we have demonstrated that a QSE-reduced a network can be 
used as a replacement for full 14 element a network when modeling silicon 
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burning, without significant errors in energy generation or nucleosynthesis. 
For such a small system, the computational benefits of such equilibrium net- 
work hybrids are also small, a factor of 2 in network computational time (due 
in most part to the factor of 2 reduction in matrix size) and a factor of 2 
reduction in the number of nuclear variables which must be evolved within 
a hydrodynamic model. For larger nuclear networks, like the realistic hybrid 
networks for the r-process [36], rp-process [47] and silicon burning [136], the 
potential for improvement in speed and size is even greater. In comparison to 
full networks with hundreds of nuclei, reduction in the number of indepen- 
dent nuclei by a factor of 2-4 can result in increases in network speed of a 
factor of 5-10 because of the nonlinear relation between matrix size and the 
length of time to solve a matrix equation. Equally important is the reduc- 
tion in the number of nuclear abundances that are hydrodynamically evolved, 
significantly reducing the memory footprint of such calculations. 



7 Conclusion 



While it is true that our basic picture of the origin of the elements seems well 
validated, a great deal of refinement and investigation remains. Due to the 
extreme conditions required for thermonuclear reactions, nuclear processes in 
astrophysics most frequently occur in deep gravitational potential wells, often 
strongly obscured from view. As a result, the number of direct observables 
are few, and nucleosynthetic observables depend strongly on the evolution of 
the stars and stellar explosions which produce them. As each of these related 
fields becomes more sophisticated, greater demands of speed and accuracy are 
placed on nuclear astrophysics calculations, mandating continued improve- 
ments in the numerical methods used. The best example of this is the recent 
trend toward mult i- dimensional hydrodynamic simulations, which greatly in- 
creases the number of independent nucleosynthesis calculations which must 
be performed. 

In this article, we have concentrated on the numerical methods of use for nucle- 
osynthesis, describing the important means by which nuclear compositions are 
currently evolved within an astrophysical simulation, thermonuclear rate equa- 
tions and Nuclear Statistical Equilibrium. While NSE solutions are much more 
economical, principally because of the much smaller number of free parame- 
ters which must be evolved, they are applicable to only a few situations. The 
principle difficulty with evolution via rate equation is the computational cost, 
which results primarily from the stiffness of the system of nuclear reactions. 
This extreme stiffness requires implicit solution and has thus far generally 
precluded the use of integration methods which do not rely on the Jacobian 
matrix. Such non- Jacobian methods remain highly sought after as a means to 
reduce the computational cost of nucleosynthesis calculations. For Jacobian 
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based integration methods, there remain considerable economies to be gained 
by taking better advantage of the sparse nature of the Jacobian. 

Physically motivated approaches can also be extremely useful in reducing the 
computational cost. One such is the use of local equilibria to reduce the size of 
the system of rate equations (and reduce problems with matrix singularity). 
Methods based on local equilibria are applicable to many situations where 
global equilibrium has not been achieved. Though the use of hybrid equilib- 
rium networks is in its infancy, it seems that in most of the situations where 
one would have heretofore used a large network coupled with hydrodynamics, 
the hybrid equilibrium networks provide sufficient accuracy and considerable 
reduction in the computational cost. 
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